##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(ggplot2)
library(ggsci)
library(ggrepel)
library(ggpubr)
library(parallel)
library(limma)

##########################################################################################
option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--input_file"), type = "character"),
    make_option(c("--width"), type = "character"),
    make_option(c("--height"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ## 突变信号文件
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"

    ## 免疫浸润
    input_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/hallmarks_gsva/gsva.tsv"    

    ## 输出
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/hallmarks_gsva"

    width <- width
    height <- height

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
input_file <- opt$input_file
out_path <- opt$out_path
height <- as.numeric(opt$height)
width <- as.numeric(opt$width)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_cibersort <- data.frame(fread(input_file))

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

rownames(dat_cibersort) <- dat_cibersort$Sample
sampleName <- rownames(dat_cibersort)
pathName <- colnames(dat_cibersort)[2]

if(ncol(dat_cibersort) > 2){
    dat_cibersort <- t(dat_cibersort[,2:ncol(dat_cibersort)])
}else{
    dat_cibersort <- data.frame(t(dat_cibersort[,2]))
    colnames(dat_cibersort) <- sampleName
    rownames(dat_cibersort) <- pathName
}

##########################################################################################
## 按人的不同病理类型合并样本
## 若一个人同一病理类型多个样本，均中位数
dat_tpm_all <- Reduce(function(x,y)bind_cols( x , y),mclapply(unique(info$ID) , function(id){

    tmp_info <- subset( info , ID == id )

    result <- data.frame()
    ## 若一个人同一病理类型多个样本，均中位数
    for(class in unique(tmp_info$Class)){
        tmp_sample <- tmp_info[tmp_info$Class==class,"Tumor"]
        tmp_tpm <- dat_cibersort[,tmp_sample]

        if(length(tmp_tpm) > nrow(dat_cibersort)){
            value <- apply( tmp_tpm , 1 , median )
        }else{
            value <- tmp_tpm
        }
        
        result_tmp <- data.frame( sample = value )
        colnames(result_tmp)[1] <- paste0( id , "_" , class)

        if( nrow(result) > 0){
            result <- cbind(result , result_tmp)
        }else{
            result <- result_tmp
        }
    }

    ## Normal样本
    tmp_sample <- unique(tmp_info$Normal)
    if(tmp_sample!="#N/A"){
        tmp_tpm <- dat_cibersort[,tmp_sample]
        value <- tmp_tpm
        result_tmp <- data.frame( sample = value )
        colnames(result_tmp)[1] <- paste0( id , "_Normal")

        ## 输出结果
        result <- cbind(result , result_tmp)
    }
    
    result

},mc.cores=10))

## 只关注Normal、IM、IGC、DGC的表达情况
col_names <- grep( paste( class_type , collapse="|") , colnames(dat_tpm_all) , value = T )
dat_tpm_all <- dat_tpm_all[,col_names]

##########################################################################################
## 与Normal做差异，使用limma包
## gsva有正负，代表了基因集合相比于不在里面的基因集，正向还是反向富集
## https://cloud.tencent.com/developer/article/1729344
## https://www.jianshu.com/p/1a82b63d65c3
result_diff <- c()
for( j in 2:length(class_type) ){
    
    class1 <- "Normal"
    class2 <- class_type[j]

    class1_col <- grep( class1 , colnames(dat_tpm_all) )
    class2_col <- grep( class2 , colnames(dat_tpm_all) )
    es_max <- dat_tpm_all[,c(class1_col , class2_col)]

    ## 保证foldchange是class2 vs class1
    group_list <- factor(sapply(strsplit(colnames(es_max) , "_") , "[" , 2) , levels = c(class1 , class2) , order = T) 
    design <- model.matrix(~0+factor(group_list))
    colnames(design)=levels(factor(group_list))
    rownames(design)=colnames(es_max)
    contrast.matrix <- makeContrasts(paste0(c(class2 , class1),collapse = "-"),levels = design)

    fit <- lmFit(es_max,design)
    ##step2
    fit2 <- contrasts.fit(fit, contrast.matrix) 
    fit2 <- eBayes(fit2)  ## default no trend !!!
    tempOutput = topTable(fit2, coef=1, n=Inf)

    tempOutput$pathName <- rownames(tempOutput)

    ## class1代表分子
    tempOutput$class1 <- class2
    tempOutput$class2 <- class1

    result_diff <- rbind( result_diff , tempOutput)
}


image_name <- paste0( out_path , "/gsva.plot.tsv" )

write.table( result_diff , image_name , row.names = F , quote = F , sep = "\t"  )

##########################################################################################
## 显著通路
dat_plot <- result_diff
dat_plot$pathway <- gsub( "HALLMARK_" , "" , dat_plot$pathName )
dat_plot$pathway <- gsub( "KEGG_" , "" , dat_plot$pathway )
dat_plot$pathway <- gsub( "MANUAL_" , "" , dat_plot$pathway )
dat_plot$pathway <- gsub( "GOBP_" , "" , dat_plot$pathway )

dat_plot$class1 <- factor( dat_plot$class1 , levels = c("IM" , "IGC" , "DGC") )
dat_plot$direction <- "nochange"
dat_plot$direction <- ifelse(dat_plot$logFC > 0 & dat_plot$adj.P.Val < 0.05  , "Up" , dat_plot$direction )
dat_plot$direction <- ifelse(dat_plot$logFC < 0 & dat_plot$adj.P.Val < 0.05  , "Down" , dat_plot$direction )
dat_plot$abs_logfc <- 2 ^ abs(dat_plot$logFC)
pathway_order <- unique(dat_plot[order(dat_plot$abs_logfc),"pathway"])
dat_plot$pathway <- factor( dat_plot$pathway , levels = pathway_order , order = T )

image_name <- paste0( out_path , "/gsva.pdf" )

col_direction <- c("blue" , "grey" , "red")
names(col_direction) <- c("Down" , "nochange" , "Up")

p <- ggplot(dat_plot, aes(x=class1, y=pathway)) + 
  geom_point(aes(size=abs_logfc,color=direction))+
  scale_color_manual(values = col_direction)+ 
  theme_bw()+  #设置背景
  labs(size="FC vs Normal") +
  xlab("")+
  ylab("")+
  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
    legend.position ='right',
    panel.grid.major=element_line(colour=NA),
    plot.title = element_text(size = 8,color="black",face='bold'),
    legend.text = element_text(size = 10,color="black",face='bold'),
    axis.text.y = element_text(size = 10,color="black",face='bold'),
    axis.title.x = element_text(size = 10,color="black",face='bold'),
    axis.title.y = element_text(size = 10,color="black",face='bold'),
    axis.ticks.x = element_blank(),
    axis.text.x = element_text(size = 10,color="black",face='bold') ,
    axis.line = element_line(size = 0.5)) 
ggsave( image_name , p , height=height, width=width) 
